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Abstract 

We consider the wave equation with highly oscillatory initial data, where there is 
uncertainty in the wave speed, initial phase and/or initial amplitude. To estimate quan¬ 
tities of interest related to the solution and their statistics, we combine a high-frequency 
method based on Gaussian beams with sparse stochastic collocation. Although the 
wave solution, is highly oscillatory in both physical and stochastic spaces, we pro¬ 
vide theoretical arguments and numerical evidence that quantities of interest based on 
local averages of \u^\‘^ are smooth, with derivatives in the stochastic space uniformly 
bounded in s, where s denotes the short wavelength. This observable related regular¬ 
ity makes the sparse stochastic collocation approach more efficient than Monte Carlo 
methods. We present numerical tests that demonstrate this advantage. 

1 Introduction 

The propagation of high-frequency waves is modeled by hyperbolic partial differential equa¬ 
tions (PDEs) with highly oscillatory solutions. In these models, the wavelength 6: is very short 
compared to variations in the structure of the medium and the wave propagation distance, 
resulting in multiscale problems. In addition, these models are often subject to uncertainty 
due to both incomplete knowledge of the model’s parameters (such as wave speed) and the 
intrinsic variability of the physical system (such as the location of an earthquake’s hypocen- 
ter or the time frequency of volcanic forces). The problem therefore has two components: 
multiple scales and uncertainty. High levels of confidence in the predictions require under¬ 
standing of the uncertainties in the model. This understanding can be obtained by a process 
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called uncertainty quantification (UQ). Furthermore, the systematic coupling and interaction 
of multiple scales and uncertainty must be considered. 

Despite recent advances in the uncertainty quantification of a wide range of mathematical 
models, see, e.g., [iii[inii2sii2ii[sn], viable UQ methodologies for high-frequency waves are 
not well developed. Current methods are based on Monte Carlo sampling techniques [2H], 
which feature slow convergence rates. The main reason that more efficient methods with 
fast spectral convergence rates are not developed for high-frequency waves is that the wave 
solution is highly oscillatory and consequently its derivatives in the stochastic space are not 
bounded uniformly with respect to £. Spectral techniques in the stochastic space have to 
resolve these oscillations and will thus require very many sample points to be accurate. Here, 
we aim at developing efficient computational techniques that feature fast convergence rates 
independent of the oscillations. 

In this work, we are concerned mainly with the forward propagation of uncertainty in 
high-frequency waves, where the uncertainty in the input parameters, such as wave speed 
and initial data, propagates through the multiscale hyperbolic model to give information 
about uncertain output quantities of interest (Qols). As the prototype model equation for 
high-frequency waves, we consider the scalar wave equation with highly oscillatory initial 
data. The main sources of uncertainty are the wave speed and/or the initial data, which are 
here described by a finite number of independent random variables with known probability 
distributions. We propose a novel stochastic spectral asymptotic method, which combines 
two techniques: (1) a Gaussian beam summation method for propagating high frequency 
waves; and (2) a sparse stochastic collocation method for propagating uncertainty and ap¬ 
proximating the statistics of Qols. An important property is the stochastic regularity of the 
Qols. By this, we mean the regularity of output Qols with respect to input random pa¬ 
rameters. The fast spectral convergence of the proposed method depends crucially on the 
presence of high stochastic regularity, independent of the wave frequency s. In particular, 
through both theoretical arguments for simplified problems and numerical experiments for 
more complicated problems, we show that although the derivatives of highly oscillatory wave 
solution u® with respect to the random parameters are not bounded independently of e, phys¬ 
ically motivated Qols based on local averages of are smooth with uniformly bounded 
derivatives in the stochastic space. Consequently, despite the slow algebraic convergence of 
the wave solution, the proposed method gives fast spectral convergence for those Qols. 

We note that the two methods composing the proposed method have already been em¬ 
ployed and are not new. The stochastic collocation method has been widely used in forward 
propagation of uncertainty in many PDF models, see, e.g., 

The Gaussian beam summation method has also been successfully applied to deterministic 
high-frequency wave propagation problems, see, e.g., [221 [23 123 123 113 [13 EZ]- However, 
combining the two methods is not straightforward since the convergence rate of the resulting 
method strongly depends on the stochastic regularity of the Qol, which in general depends 
on £ and in principle may be as bad as the regularity of u®. The main contributions of this 
work include: (1) constructing a novel algorithm for the uncertainty propagation of high- 
frequency waves by systematic coupling of sparse stochastic collocation and Gaussian beam 
approximation; and (2) verifying the fast spectral convergence of the proposed algorithm for 
particular types of quadratic quantities, independent of frequency. 
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The rest of the work is organized as follows. In Section 2, we formulate the high-frequency 
wave propagation problem with stochastic parameters. The stochastic spectral asymptotic 
method for the forward propagation of uncertainty in the problem is presented in Section 
3. In Section 4, we discuss the stochastic regularity of the highly oscillatory solution and 
a particular quantity of interest, through theoretical arguments and numerical experiments. 
In Section 5, we provide numerical examples and demonstrate the spectral convergence of 
the proposed method for the quantity of interest. Finally, we summarize our conclusions in 
Section 6. 

2 Problem Statement 

Consider the Cauchy problem for the stochastic wave equation 

«tt(Cx,y)-c(x,y)2 AM^(t,x,y) = 0, in [0, T] x E** x T, (la) 

«®(0,x,y) = 5 (^(x,y), uf(0, x, y) = 5 f^(x, y), on {t = 0} x M” x T, (lb) 

where e C is the stochastic wave function, t e [0,T] is the time, x = (xi,... ,x„) e E”' 
is the spatial variable, y = (r/i,... ^Vn) G T C E"^ is a random vector, and £ <C 1 is a 
small positive number, representing a typical short wavelength defined below. We use the 
convention that V represents the gradient with respect to the spatial variables x. 

Sources of uncertainty. The uncertainty in model 0 is described by a random vector, 
y, consisting of A e N+ independent random variables, yi, ... ,?/Ar, with a bounded joint 
probability density, p(y) = PniVn) ■ T —)■ E+. There may be two sources of uncertainty: 
uncertainty in the local wave speed, c = c(x, y), and uncertainty in the initial data, gf = 
pf(x,y) and = 5 '|(x, y). We make the following smoothness, uniform coercivity and 
boundedness assumptions on the wave speed: 

ci-,y)eC^iR^), VyGT, (2) 

0 < Cmin < c(x,y) < Cmax < OO, VxgE”, VyGF. (3) 

In the case when T is not compact (for example when random variables are Gaussian), in 
addition to ([^ , we assume the spatial derivatives of the wave speed are uniformly bounded in 
r. ifF is compact (for example when random variables are uniform), this extra assumption is 
automatically satisfied. We further assume that the random wave speed has bounded mixed 
y-derivatives of any order, i.e., for a multi-index k= {ki,... , fcjv) G with |k| > 0, we 
define := k., and assume 

... Uyjyj 


dyc{ •, y) 




< oo, VyGF. 


High-frequency waves. We consider highly oscillatory initial data of the form 


Pf(x,y) = ^(x,y) e 


,i $o(x,y)/£ 


^2(x,y) = j5o(x,y) 


i$o(x,y)/£ 


(4) 

(5) 
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where the short wavelength, £, is assumed to be much smaller than the typical scale of the 
structure of the medium (variations in the wave speed) and the wave propagation distance 
(the size of the computational domain). Such initial data generate high-frequency waves prop¬ 
agating in low-frequency media. The functions Aq, Bq and $0 are assumed to be real-valued 
and smooth, with |V$| bounded away from zero, Vy G T. Consequently, we consistently 
assume 

P 2 (->y) e ^^(M”), VyGT, V£>0. (6) 

Well-posedness of the problem. By formally partial-differentiating problem 0 with 
respect to t and x and a simple extension of the proof for deterministic problems nziiiaEi, 
we can show that problem 0 is well-posed. In other words, with the random wave speed 
satisfying 0-0 and the initial data satisfying ([^, there exists a unique solution u^( • ,y) G 
C~([0,T] X M”) for every y G T to problem Q, which depends continuously on the data. 
Moreover, at a fixed (t, x), the solution belongs to the Hilbert space of stochastic functions 
with bounded second moments, i.e. 

u%t, X, .) G Lj(r) := |u : r ^ E, ^ k(y) ^ p(y) dy < ooj , 

where the space is endowed with its standard inner product. We also note that by 
partial-differentiating the problem ([^ with respect to y and employing assumption 0 , we 
can further show u® G (^^“([O, T] x E” x T). We refer to [SUES] for a more detailed discussion 
on the well-posedness of stochastic wave propagation problems. 

The ultimate goal is the prediction of statistical moments of the high oscillatory solution, 
u^, or the statistics of some given physically motivated Qols. In particular, we consider the 
following Qol 

Q^{y)=[ |u^(T,x,y)|V(x)dx, (7) 

where 'tp is a given real-valued function, referred to as the Qol test function. Throughout this 
paper, the Qol test functions we consider will always be smooth and compactly supported. 
Ip G ^^(E”), unless otherwise stated. 

Remark 1. In the present work, we concentrate only on the quadratic Qol 0. which rep¬ 
resents the average strength of the wave inside the support of pj. Other types of Qols will 
be considered and studied in future works. We note in particular that the linear quantity 
Q^{y) = /ru u^{T, X, y) 'ippx) dx typically vanishes as £ ^ 0 and is therefore not of interest. 

3 A Stochastic Spectral Asymptotic Method 

In this section, we present an efficient numerical method for solving the stochastic wave equa¬ 
tion ([^ with highly oscillatory initial functions ([^. The method combines two techniques: 
a Gaussian beam summation method for propagating high-frequency waves; and a sparse 
stochastic collocation method for approximating the statistics of Qols. 


4 


3.1 Gaussian beam approximation 

The Gaussian beam method describes high-frequency waves in a way closely related to geo¬ 
metrical optics and ray tracing. Geometrical optics formally assumes the solution of 0 to 
be of the form 

= ( 8 ) 

As the phase (f) and amplitude o are independent of frequency and vary on a much coarser 
scale than does the solution, they might be resolved with a computational cost independent 
of £. In the limit e ^ oo, we arrive at the eikonal and transport equations [33]: 


+ c(x,y)|V(/)| =0, 


Qjf 


c2(x,y)A0 - c(x,y)Va • 


2c(x,y)|V(/)| |V0| 

where (t, x, y) G (0, T] x E” x T. Initial data is given as 


= 0 , 


(9) 

( 10 ) 


a(0, X, y) = Ao(x, y), 0(0, x, y) = d>o(x, y), 

with Aq and in 0. 

The bicharacteristics (q(t, y), p(t, y)) of the eikonal equation a satisfy 


( 11 ) 


Here the vector q is usually called the ray and the vector p is called the slowness. For 
initial values q(0,y) = zo(y), p(0,y) = V'I>o(zo(y),y) we have that |0(t,q(t,y),y) = 
0 and thus the phase 0 is given by the simple expression 0(t, q(t, y), y) = d>o(zo(y), y). 
Moreover, the slowness vector always points in the direction of the phase gradient, p{t, y) = 
V0(t, q(t, y), y). In this way, the eikonal equation always admits a smooth solution locally 
in time, more precisely up to the point where two rays cross. If we denote the ray starting 
at z by q(t, y; z), then caustic points are those where this ray function degenerates in z, i.e., 
where det 5q(t, y; z)/9z = 0. At these points, the initial ray position fails to be a well-defined 
function of the current ray position. A global smooth solution of (9) is in general not feasible 
as caustics develop when rays cross. The form of the solution in (8) is thus not correct. The 
solution becomes multi-phased and geometrical optics predicts an unbounded amplitude [3S| , 
cf. the discussion in Section IH 

Gaussian beams present another type of high frequency approximation closely related 
to geometrical optics evaluated along rays. However, unlike geometrical optics, a Gaussian 
beam is well-defined for every t and their superposition performs well even at caustics. 

A Gaussian beam, u®(t,x,y), has the same form as (|^: 

( 12 ) 




where the phase, •F, and amplitude. A, are centered around the geometrical optics ray, q(t, y): 


X, y) = a{t, X - q(t, y), y), 
X, y) = 0(t, X - q(t, y), y). 
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The two methods differ in their assumption on the phase, <p. Geometrical optics assumes 
(p is real whereas the Gaussian beam method sets the phase to be complex away from the 
ray. The imaginary part is chosen such that the Gaussian beam (12) decreases exponentially 
away from the central ray. 

For the simplest (first order) Gaussian beams, o and (p read 


o(Tx,y) = oo(t,y), 

0(Tx,y) = (/)o(t,y) + x • p(t,y) + • M(t,y)x, 


(13) 

(14) 


where M is a symmetric matrix with a positive imaginary part. Note that these expressions 
can be interpreted as a zeroth and a second order Taylor expansion in x of the amplitude and 
the phase respectively; the slowness p is then the phase gradient, as in geometrical optics, 
while the M matrix is the Hessian of the phase. By choosing q, p, M and oq in the right 
way, the phase will satisfy the eikonal equation ([^ up to 0(|x — q(t,y)P) and the amplitude 
will satisfy the transport equation (10) up to 0(|x — q(t,y)|). The coefficients (pQ, M, ao then 
obey the following ordinary differential equations (ODE) (see [32]): 


(po — 0, 

M ^D + B^M + MB + MGM, 
1 


ao 


2|p| 


c(q, y) Tr(M) - Vc(q, y) • p - 


c(q,y) p • Mp' 


ao. 


(15a) 

(15b) 

(15c) 


with 


D = |p|W(q,y), 


B = 


p® Vc(q,y) 


^ ^ c(q,y) j c(q,y) 


We will always assume that the solution to these ODEs exists for the time intervals considered. 
A sufficient condition is, for instance, that c(x, y) and all its x-derivatives up to order three 
are uniformly bounded for all x G M” and y G F. 

Remark 2. Geometrical optics also allows for an alternative to the eikonal and transport 
equation pair and ([l0|), namely: 


(pt - c(x,y)|V(/)| = 0, 


c^(x,y)A(p - (ptt c(x,y)Va-V0_ 

Gj-t ^ / \ I _ , I G I _ , I — U • 


2c(x, y)|V,;)| 


IV./.I 


This alternative corresponds to waves moving in the opposite direction and leads to opposite 
signs in the ODEs and ( [T^ . 

The name of the Gaussian beam is substantiated by its Gaussian shape. The imaginary 
part of (p is determined by x • Im(M) x; thus 

X, y)| = ao exp - q(^, y)) • Im(M)(x - q(t, y))^. 
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Since M has a positive imaginary part, \v{t, x, y) | is Gaussian with a width |x — q(t, y) | ~ ^/e. 
The positivity of Im(M) is key to the Gaussian beam approximation. Indeed, |32] states that 
M = and Im(M) > 0 hold true at any time, provided they are valid for the initial data. 
In general, to construct higher order Gaussian beams, higher order terms in the Taylor 


expansion of 0 and a in (13), and (14) respectively, must be employed. Analogously to the 
above derivation, x) is required to solve the eikonal equation (|^ to 0(|x — q(t, y)!^’*'^) 
for the fc-th order beams and the amplitude terms aj to solve the transport equations to 
0(|x—q(t, y)!^”^-’), respectively. This again translates to a system of ODEs for the coeiircients 
in the Taylor expansions. 


To approximate more general solutions, we use a superposition of Gaussian beams 

1 




(16) 


over a compact domain, Kq C M”, where the initial data is supported, and x, y;z) is a 
Gaussian beam that starts in the point z & Kq. 


Due to linearity of the wave equation the superposition in (16) is also an asymptotic 


solution, since each individual beam is such a solution. The pre-factor (27re)~^^‘^ normalizes 
the superposition such that = 0(1) away from caustics. The coefficients of the beams 
are similarly parameterized by z and denoted q(t,y;z), p(t,y;z), etc. Following 
choose the initial values to be 


we 


q(0,y;z) = z, 

(17a) 

p(0,y;z) = V4>o(z,y), 

(17b) 

M(0, y; z) = V^$o(z, y) + *1, 

(17c) 

0(0, y;z) = $o(z,y), 

(17d) 

ao(0,y;z) = Ao(z,y). 

(17e) 


It has been shown in na that the initial data for fc-th order beams can be chosen such that 
the error in an £-scaled energy norm satisfies the estimate 

sup •) - “^(^1 OIU < C{T)e^, 

t&[0,T] 


for some constant C{T) independent of £. Note that this holds regardless of the presence of 
on. Numeric 

“Gs(^>x,y) 


caustics in the solution. Numerically, integral (16) is approximated by the trapezoidal rule: 

1 




x,y;z,) As”, 


(18) 


where Az ~ ^/£, and the ODEs ( |TT| ) and (15) are solved with a numerical ODE method. 
The computational cost of the Gaussian beam method is then much smaller than that of a 
direct solver. 
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3.2 Sparse Stochastic Collocation 

For stochastic PDF models, such as the wave model 0 with stochastic inputs, there are 
in general two types of methods for propagating uncertainty: intrusive and non-intrusive. 
Intrusive methods, such as perturbation expansion and stochastic Galerkin [m sn [221 El 
ES], require extensive modihcations of existing deterministic solvers. On the contrary, non- 
intrusive methods, such as Monte Carlo da and stochastic collocation iniEiEii, are sample- 
based approaches. They rely on a set of deterministic models corresponding to a set of 
realizations and hence require no modihcations of existing deterministic solvers. In this 
work, we consider the stochastic collocation method on sparse grids and briehy review the 
method for the uncertainty propagation of high-frequency waves. 

The stochastic collocation method consists of three main steps. First, the problem Q 
is discretized in the physical space, i.e., in space and time, using a deterministic asymptotic 
solver such as the Gaussian beam method described in Section EH We therefore obtain a 
semi-discrete solution UQB(t,x,y), keeping the variable y in the stochastic space continu¬ 
ous. The semi-discrete problem is then collocated on a set of 77 G N+ collocation points, 
^ r, to compute rj approximate solutions, UQB(t, x, Finally, a global polyno¬ 

mial approximation, Ugb,??) built upon those evaluations, 

“GB,„(^>x,y) = (19) 

k=i 

for suitable multivariate polynomials, {Tfc(y)}^^^, such as Lagrange polynomials, in the 
stochastic space. 

A key point in the stochastic collocation method is the choice of the set of collocation 
points, {y^^^}fc=i, i.e., the type of computational grid in the A^-dimensional stochastic space. 
A full tensor grid, based on the Gartesian product of mono-dimensional grids, can only be 
used when the dimension of the stochastic space, A’, is small, since the computational cost 
grows exponentially fast with N {the curse of dimensionality). To clarify this, let i G Z_|_ be 
a non-negative integer, called the level. Moreover, for a given index, j G let p{j) be an 
increasing function, which relates index j to the polynomial degree and hence to the number 
of interpolation nodes. Typical functions for the polynomial degree are given by the linear 
rule 

P{j) = j, (20) 

and the nested rule 

p{j) = 2"' for j > 0, p(0) = 0. (21) 

In the full tensor grid, we take all polynomials of degree at most p{£) in each direction, and 
T] = {p{i) + 1)"^ grid points (nodes or knots) are therefore needed. 

Alternatively, sparse grids can reduce the curse of dimensionality. They were originally 
introduced by Smolyak for high-dimensional quadrature and interpolation computations [31] . 
To understand the general sparse grid construction, let j = {ji,. ■ ■, Jn) G be a multi¬ 
index containing non-negative integers. In each direction, F^^, with n = 1,... ,N and for a 
non-negative index jn in j, we introduce the univariate polynomial interpolation operator 

W ": G»(r,.) ^ p,o.)(r„), 
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on p{jn) + 1 suitable knots. Here, Pp is the space of univariate polynomials of degree p. 
The univariate interpolation operator applied to a stochastic function of only one random 
variable, say v(yn), reads 


P(in)+1 

U^’^[v{yn)] := ^{yn^)Lk{yn)- 

k=l 

Here, Lkiyjn) is the univariate Lagrange polynomial of degree k — 1. With U~^ = 0, we then 
dehne the detail operator as 

After the introduction of a sequence of index sets, X(£) C the sparse grid approximation 
of the solution to 0 at level I reads 


N 

w®(-,y) ~‘5i(q[u^B(-,y)] = ^ [u^B(-,y)]. (22) 

iex{i) n=l 


Equivalently, we can rewrite the sparse approximation (22) as 


N 


^x{e)[uhBi • >y)] = E <^(j) [“gb( • >y)]> 

jex(e) n=l 


c(j)= y (- 1 )'" 


(23) 


where the multivariate interpolation operator reads 


N 




n=l 


P(ji)+1 P(jw)+1 

•W)]:= ••• £ «GB 

ki=l kN=l 


{■,{yi'\---,yN''^))Lk{y), 


i={0,l}^ 

i+jex(r) 




(24) 


and Tfc(y) = Yln=i ^kniyn)- This operator is given by the tensor product of N univariate 
interpolation operators, constructed on a full tensor grid corresponding to the multi-index j. 
The second formulation (23) shows that a sparse grid is a linear combination of a few tensor 
product grids, each with a relatively small number of grid points. 


To characterize the sparse approximation operator in (22) and (23) fully, we need to 
provide the following: 


(1) A level i G Z+ and a function p(j) given by either (20) or (21). 

(2) A sequence of sets X(£). Typical examples of index sets include: 

o Total degree index set: Xtd(^) = {j : Y^n=i3n < 
o Hyperbolic cross index set: Xhc(^) = {j : I{n=i{jn + 1) < ^ + 1}. 

(3) The family of grid points to be used. Typical choices include: 

o Gauss abscissas, that are the zeros of p-orthogonal polynomials; see e.g. EH. 
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o Clenshaw-Curtis abscissas, that are the extrema of Chebyshev polynomials; see 

e.g. [H7]. 

o Leja abscissas; see e.g. j2H]- 

Remark 3. (Full tensor and Smolyak grids) The full tensor grid eorresponds to the full 
tensor index set, = {j : max„ < F\. The Smolyak sparse grid is a particular type of 

sparse grid, where we use the nested rule ( plj ) and the total degree index set, based on either 
Gauss or Clenshaw-Curtis abscissas. In the latter case, we obtain a nested grid, in which 
grids in different levels are embedded. 

In many practical applications, the main goal is the computation of the statistical mo¬ 
ments of the solution or some Qols. In such cases, we can directly compute the statistical 
moments using Gauss or Clenshaw-Curtis quadrature formulas, without explicitly construct¬ 
ing the solution or other quantities by sparse approximation formulas such as ( |^ . Suppose 
we want to compute the statistics of an operator applied to the solution, say F{u^). This can 
for instance be the solution itself, the Qol Q, or different powers of these quantities when 
computing higher moments. We write 


E[F(u®(-,y))] RiE[5i(q[F(u^B(-,y))]] = / Sx{i)[F{uf^{-,y))]p{y) dy 

k=i 


where, by (23)-(24), the quadrature weights read 

Ok = Ck j^Lkfy) p{y)dy. 


Here, each global index, k = 1,... ,p, corresponds to a local multi-index [ki,..., k^] in the 
formulas (p3|)-p4|). Moreover, the coefficients Ck correspond to the coefficients C in (1^. 


Convergence of stochastic collocation. It is well known that the rate of convergence of 
stochastic collocation for a stochastic function, say v(y), in general depends on the stochastic 
regularity of the function, i.e., the regularity of the mapping u : T —)■ M. Fast convergence 
is attained in the presence of high stochastic regularity. For instance, suppose that v{y) is 
continuous and admits an analytic extension in the complex region 


fr = = {zi,..., Z]\r) G I for some j: dist(2;j, Fj) < r and .2^ G F^ for fc 7^ , 

whose size is characterized by the radius of analyticity, r > 0. Then, the maximum error, 
Cmax; in the sparse approximation, <Sx(q[u(y)], in (23) on the Smolyak grid based on Clenshaw- 
Curtis abcissas satisfies (see Theorem 3.10 in BB) 


fmax <C(T.N.r) M{V. T, T) J/ 


•f7/(l+log(2 AT)) 

5 


M{v, T, F) = max |n(z)|, 

zGF T 


where a > 0 is directly proportional to the radius of analyticity, r, and C, M are independent 
of rj. Therefore, the larger the radius of analyticity of the map, the faster the convergence rate 
in rj; see also In the case of high-frequency waves with highly oscillatory solutions, 

r, a and M in general all depend on the wavelength 6:. To maintain fast convergence, it 
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is therefore important that M and 1/r remain bounded as £ ^ 0. Since the size of M 
depends on the size of the y-derivatives of u on F, we need uniform bounds for y-derivatives 
independent of £ to obtain fast convergence for all £. In other words, we require high regularity 
of the mapping : F ^ M, uniformly in £. This poses extra challenges in employing 
stochastic spectral techniques for high-frequency waves. The stochastic regularity and the 
convergence of stochastic collocation for high frequency waves is discussed in more detail in 
Section IH 


4 Stochastic Regularity of High-Frequency Waves 

For fast convergence, regular dependence of the quantity of interest on the input random 
variables is required, as discussed in Section 3.2, and, ideally, Q^(y) G C'°°(F). For high- 


frequency problems, we would like to have an even stronger bound, namely 


sup 

yer 


d'^Q^iy) 


dy^ 


< Ck, Vk G 


+ ’ 


(26) 


where the constants (Fk can be taken independent of the wavelength, £. If this is not true, 
the spectral convergence rate of a stochastic collocation method could deteriorate for small 
wavelengths, £; with the bound (25), we expect that there is a rate that is uniform in £. 

We note that in general M®(t,x,y) will be oscillating with period ~ £ in both x and y. 


If the corresponding quantity of interest would inherit this property, a bound like (25) could 


not hold. Nevertheless, in this section, we conjecture that for a Gaussian beam superposition 
approximation of with the initial data given by ([^ and (15) respectively, the estimate 
(iMl) holds. 


Conjecture 1. Bound (25) holds for the following quantity of interest computed with the 
Gaussian beam approximation 


Q’(j)=f |MGB(r,x,y)|V(x)<(x, 

JR" 


ifi^e 


This is based on formal theoretical arguments and several numerical experiments, which 
are presented in the subsections below. The conjecture will be proved in an upcoming work. 


4.1 Motivation away from caustics 


In this section, we give a non-rigorous argument for why (25) should hold, at least in the 


case when there are no caustics in the support of the quantity of interest. We thus consider 
a high-frequency solution for 0 and we assume (to, xq, yo) G [0, T] x M”' x F is a non-caustic 
point. By the theory of Maslov 1201E] and Fourier Integral Operators [nmaiz] there is 
then at least a neighborhood, [/, of (to,Xo,yo) in which the solution has the form 


K 


u 


(t,x,y) = + V(i.x.y) e U, 


(26) 


i=i 
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for some integers K, which represents the number of waves passing through (to,xo), and mj, 
which is the Keller-Maslov index of wave j, i.e., the number of times it has passed through a 
caustic. Note that in general K > 1 when to > 0, even if the initial data is of the single wave 
form ([^ since new terms appear when caustics develop; see ra for further explanations. 
The amplitudes Aj{t,:K,y) and phases 0j(t,x,y) satisfy the geometrical optics equations, 
i.e., the transport and eikonal equations Q, and (10) respectively. They are independent 
of e and smooth in all variables, given the assumptions of a smooth wave speed and initial 
data. Moreover, the phases, are related to the initial phase, d>o, and to p(t,y;z) via the 
ray function, q(t,y;z), as follows. For each (t, x, y) G U and each 1 < j < K, there is a 
Zj G Zj{t,y) such that 


x = q(t,y;zj), x, y) = $o(zj, y), x, y) = p(t, y; z^), 


(27) 


where Zj{t,y) are disjoint open subsets of supp(Ao( • ,y))- 

Based on this, let us assume that in the support of the Qol test function, xp, there are 
no caustic points at t = T, and that the solution satisfies (26) and (27). The quantity of 
interest 0 can then be written 


Q%y)=j |u"(r,x,y)|V(x)dx 


/ 

Jw 


K 


Y, Aj (T, X, y)e*'^i(Tx,y)/e+im,-^/4 

K K 


'0(x)(ix 




i=i j=i 

where 5? denotes the real part and 


Qljiy)= [ A(2",x,y)^j(r,x,y)e*(^' 






(28) 


We note that (25) holds if the same estimate holds for each Qlj{y). Clearly the diagonal 
terms, 

are smooth and independent of e, hence satisfying ( p5| ). For the off-diagonal entries, i ^ j, 
however, we have for a multi-index, k G ZN, 


/? 1 S JRP \ S I 


ayk 


* ^ J, 


e=i 


where are smooth, e-independent functions built from sums and products of y-derivatives 
of Ai, Aj, (pi and (pj evaluated at t = T. Because of the pre-factors, e“^, it follows that we 
cannot bound as in (25) unless 


Jw^ \e 


Vm G Z_i 


(29) 
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To verify this, we use stationary phase arguments. More specifically, the following theorem 
formalizes the dependence of oscillatory integrals on critical points of the phase; see [13]. 


Theorem 1 (Principle of (non)-stationary phase). Let G and ip G 

such that Vip ^ 0 on supp{ip). Then, for all m G there exist constants, Cm, independent 

of e, such that 





< Cms^. 


When applying this theorem to the integral in ( |29| ), we thus obtain the required estimate 
provided two conditions are fulfilled. First, the Qol test function should be smooth with 
compact support, i.e., if G Since the functions, are smooth, the theorem 

applies. However, we can, for instance, not expect (25) to hold if we define the Qol just as 
an integral over a domain, D CMT, 


Q%y) = [ l«®(7',x,y)p(ix, 

JD 


(30) 


which corresponds to choosing ip as the (non-smooth) characteristic function of D. Second, 
there should be no stationary points in the support of ip. However, under our assumptions, 
this is not possible. Indeed, suppose there is a stationary point, i.e.. 


V(/>i(r,x,y) = V(/)j(T,x,y) 


for some x G supp(';/^). Then, by (27), there are two geometrical optics rays starting at 
Zj G Zi{T,y) and zj G Zj{T,y), respectively, such that 

q(T, y; z^) = q(T, y; Zj) = x, 

p(T,y;Zi) = V(/)*(T,x,y) = \/(pj{T,x,y) = p(T,y;Zj). 


But, by uniqueness of the solution to the Hamiltonian system (11), we must therefore have 


Zj = Zj. Hence, the two rays are the same and since the sets {Zj{t,y)} are disjoint, we get 
i = j, a. contradiction. 

In conclusion, the formal arguments above suggest that if the Qol test function is smooth, 


then, at least away from caustics, the bound (25) holds. 


4.2 Numerical justification and preliminaries 

In this section, we consider a number of problems with high-frequency solutions to 0 , to 
justify the bound (25). We recognize that the theoretical arguments for (25) in the previous 


section were only formal. In general high-frequency solutions, there are caustics, which we 
did not account for. Furthermore, we did not consider the regularity of the approximation 


error, 0{e), in (26). 

Here we will show three numerical examples to demonstrate that the quantity of interest, 
Q^, can also be smooth in more general cases. First, we show a simple one-dimensional ex¬ 
ample with constant coefficients, where the arguments in Section |4T] are essentially rigorous; 
in one dimension, there are no caustics, and for constant coefficients, geometrical optics is 
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exact. Our second example is still one-dimensional with no caustics, but the speed of prop¬ 
agation varies in space. Geometrical optics is therefore accurate only up to 0{e). The last 
example is a two-dimensional problem where there are caustics in the solution. As we will 
see, for all these cases, the Qol and its derivatives remain non-oscillatory when the Qol test 
function is smooth. 

All three examples use an initial condition of the single wave type, 

u^(0,x,y)=Ao(x,y)e'‘^°("’^)/^ (31) 

where the amplitude, Aq, is a sum of two pulses moving towards each other. The integration 
via the trapezoidal rule ( |I^ is carried out to obtain the approximated Gaussian beam solu¬ 
tion, Uqq. The spacing between the starting points is set to Az = ^/e in all examples. The 
Qol test function in 0 is based on the smooth function 

_ ( 

^(x) = I e |x| < 1, (32) 

y 0, otherwise. 

To show the regularity of at high frequencies, it is computed for a sequence of small 
wavelengths, e = (1/40,1/80,1/160). The final time in Q is fixed at r = 1 for all three 
examples below. 

The integration over variable x in 0 is again carried out by the trapezoidal rule with 
a spatial step that is uniform in all dimensions. Ax = unless stated otherwise. This 
corresponds to ten points per wavelength, which gives a very high accuracy when G (M**) 
due to the spectral convergence of the trapezoidal rule for such functions. 

4.2.1 Example 1: Constant speed of propagation in ID 

We start with the simple case of the one-dimensional wave equation 0 with constant speed 
of propagation, where we have an explicit expression for the solution u(t, x, y). We investigate 
a case when two high-frequency pulses are moving towards each other, 

Aq{x) ^ g{x - si) + g{x - S2), 5f(x) = ^o{x)=x‘^, 

where si = —S 2 = 1-5. The initial data has no uncertainty, i.e., it does not depend on y. 
However, the constant (in x) speed is uncertain, c = c{ij). The initial time-derivative of u is 
taken such that the pulses propagate towards each other, 

( 0 , X, y) = -c{y) {^'{x - sQ -F g{x - sQ j 

-6 c{y) ^g'{x - S2) + g{x - 82)^ 

Moreover, we let c{y) = y > 0. The solution is then given by d’Alembert as 

U%t, X, y) = g{x - Si - yt)e^'^o(x-yt)/e ^ ^(^ _ ^ y^y<Po(x+yt)/e_ (33) 
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Figure 1: Example 1. Absolute value of solution for various times, t, when y = 2 and 
£ = 1/40. The two Qol functions, 'tpo and 'ipi, are overlaid. 
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Figure 2: Example 1: Qols, Qq (top row) and Qf (bottom row), and their first and second 
derivatives, for different wavelengths, £. For Qf the curves are almost independent of £ while 
the non-smooth Qq exhibits a clear loss of regularity for decreasing e. 


The absolute value of a particular solution with e = 1/40 and y = 2 is shown at various times 
in Figure [T] 

We consider two different qnantities of interest. 


Qj(y) = j \u%T,x,y)\‘^i>j{x)dx, 


T=l, J = 0,1, 


where ipi{x) = ip{2x) is smooth, based on (32), and ipo{x) is the (non-smooth) characteristic 
function for the interval [—1/2,1/2], which reduces the form of Qq to (j^. The Qol test 
functions, i/jjix), are plotted together with the solution in Fignre[^ 

By nsing the explicit solution (33), we can write Qf{y) as a sum in the same way as in 


Section 4.1 


Qiiy) = MQIM + Quiy) + QIM + Q% 2 {y) , 


where Q\i{y) and Q\ 2 {y) independent of e and 

^ f . T>q(cc —r/) —T>Q(cc-|-r/) 

Quiy) = Qkiy) = j dix - Si - y) g{x - S2 + y)e^ ^ 'il;j{x)dx 

—4:xyi\ 


x — Si — y) g{x — S 2 + y) exp ^ xpj(x)dx. 
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r\j 


From this expression, one can derive that for 'tpQ, the derivatives behave as 


dyk 


Hence, condition (29) is not fulfilled and we should observe oscillatory behavior, 
other hand, the function G suppresses the oscillations when y 


> 


On the 
0 according to 

Theorem 4.1 and we should expect smooth higher derivatives of Q\. 

To demonstrate this numerically, we plot in Figure]^ the quantities of interest, Qq and 
Q\, with their first and second derivatives. Indeed, the derivatives of Q\ are smooth, while 
those of Qq oscillate; in particular, the second derivative of Qq grows as ^ when e 0, 
which means that (25) is violated. For accurate approximation of Qq, we use 50 points per 
wavelength to discretize the spatial grid x which makes the quadrature error negligible. 


4.2.2 Example 2: Variable speed of propagation in ID 


Problem ([^ is only analytically solvable in the d’Alembert form (33) if the speed of prop¬ 
agation c is constant. In more general cases, the Gaussian beam approximation has to be 
implemented. Analogously to the previous example, we choose the initial data as 


Ao{x, y) = g{x - Si(y)) + g{x - S2(y)), gix) = e 


-10a;2 


d>o(a:) = \x\, 


which represents two high-frequency pulses moving towards each other. We set the Gaussian 
beam initial data as in 0 The initial position of the pulses and the non-constant speed of 
propagation depend on two stochastic variables, y = {yi,y 2 )- More precisely, 

Si(y) = -S 2 (y) = 2 / 1 , c{x, y) = 1 + i + 2 / 26 ^"’+^^") . 

The absolute value of the solution with e = 1/80 for y = (1,0) at different times is plotted 
in Figure Figure |4] shows four different realizations of the absolute value of the solution 
with £ = 1/80 at time T = 1. In all realizations, we keep yi = 1 fixed and vary 1 / 2 • This 
corresponds to different realizations of the wave speed. The Qol test function, 'tp{x) = 'tp{2x), 
and the wave speed, c{x,y), corresponding to different realizations are also shown in the 
figure. 

In Figure we plot the quantity of interest ([^ and its first and second derivatives along 
the line y(r) = (1 -|- r, 1 — 2r). We observe that Q^ and its derivatives are smooth and do 


not oscillate with e. Hence, the estimate (25) is fulfilled. 


4.2.3 Example 3: Caustics in 2D 

In our final example, we consider a two-dimensional case with two initial wave pulses, 
^o(x,y) = p(x-si(y))+^(x-S2(y)), ^(x) = e"^W'. 

The deterministic initial phase 

4>o(x) = |a:i|-I-X 2 , ii={xi,X 2 ), 

is chosen such that a cusp caustic develops at t = 0.5, and two fold caustics form at t > 0.5. 
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Figure 3: Example 2. Absolute value of solution for t = 0.25, 0.5, 0,75, 1 when £ = 1/80 
and y = (1,0). The Qol test function, and speed, c, are overlaid. 
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Figure 4: Example 2: Four realizations of the absolute value of the solution at time T = 1 
with £ = 1/80. The Qol test function, and speed, c, are overlaid. 



(a) Q^(y(r)) 


(b) iQ^yir)) 



0.5 -“o 


(c) ^Q^(y(r)) 


Figure 5: Example 2: Qol and its hrst and second derivatives along the line y(r) = 
(1 + r, 1 — 2r) for r G [0, 0.5] and different wavelengths £. 
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(a) t = 0 



(b) t = 0.2 


(c) t = 0.4 




r 



(d) t = 0.6 


(e) t = 0.8 


(f) t = 1.0 


Figure 6: Example 3: Absolute value of solution for various times, t, when y = (1,1). 
Canstics appear for t > 0.5. Circle indicates the snpport of the Qol test fnnction. 


Both the initial location of the pulses and the constant speed of propagation are random, 
depending on two stochastic variables, y = (?/i,|/ 2 ); as 

si(y) = -S 2 (y) = (|/i,0), c(y) = y 2 . 

We consider the Qol Q with '^(x) = ip{2x.). 

In Figure|^the absolnte value of the solution at varions times is shown. In this simulation, 
the wavelength is £ = 1/40, the pulse centers are Si = —S 2 = (1,0), and the speed of 
propagation c = 1. The bold line represents the support of the Qol test function, 'tp. Figure 
shows the Qol and its derivatives along the line y(r) = (1 + r, 1 + 2r). Note that for most 
of these y-valnes, the two pnlses overlap at the final time, T — 1. 


5 Numerical Examples 


In this section, we present two nnmerical examples to demonstrate the efficiency and appli¬ 
cability of the method proposed in Section 

We consider the Canchy problem Q in a two-dimensional spatial space and let x = 
(xi,X 2 ) G We employ the proposed stochastic spectral asymptotic method to approxi¬ 
mate the solution, u®, and the expected value of the Qol in ([^. The Qol test functions are 
given in terms of the smooth function, pj G C'/°(E^), in ( |^ as 'ip{xi,X2) = pj{ 2 xi, 2 x 2 ) and 
'ip{xi,X 2 ) = ipixi — 1,X2) for the first and second numerical examples, respectively. In both 
examples, we nse the Smolyak sparse grid based on Clenshaw-Cnrtis abscissas and the nested 
rule (21). We show that fast convergence rates are obtained as predicted in Section]^ For 
each example, we also compare the convergence rate with the average rate obtained from ten 
independent Monte Carlo simulations. 

As in Section 4.2, the step size used in the quadrature approximation of 0 is chosen 


as 


Ax = ^ and the space between the Gaussian beams in (18) as Az = ^fe. 
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(a) Q"(y(r)) (b) ^Q"(y(r)) (c) |iQ^(y(r)) 


Figure 7: Example 3: Qol Q^, and its first and second derivatives along the line y(r) = 
(1 + r, 1 + 2r) where r e [0,0.5], for different wavelengths, £. 


5.1 Numerical test 1: Two pulses 


In this example, both wave speed and initial data are uncertain and described by a ran¬ 
dom vector, y = (yi,..., r/s), containing N = b independent uniformly distributed random 
variables. The constant random wave speed is given by 


c(y) = yi ^ W(0.8,1), 


and the random initial data are given by Q with 

$o(x) = |a:i|, >lo(x,y) = - si,di(y)^ + - S 2 (y),d 2 (y)^, 

where 

c/(x,d) = 

and 

Si = (-1,0), S 2 (y) = (?/2,1/3), di(y) = (2/4,5), d 2 (y) = (?/5, 2 / 5 ), 

and 

?/2-Z7(1,1.5), ?/3-Z^(0,0.5), y4-Z^(5,10), 10 ). 


Hence, the initial solntion consists of two Gaussian wave pnlses, and the vectors Sj and dj, 
with j = 1,2, represent the position and shape of the pulses, respectively. 

Figure|^shows six realizations of the magnitnde of the approximate solution \uq^{T, x, y)| 
with wavelength £ = 1/40 at the time T = 1. In each realization, the central circle indicates 
the support of the Qol test function and the other two circles/ellipses indicate the supports 
of the initial solntion that consists of two Gaussian pnlses. 

We make a convergence stndy for a set of wavelengths, £ = 1/40,1/80,1/160. For each 
wavelength, we consider different levels, t' > 1, and compnte the relative error in the expected 
value of the Qol in Q: 


mi)) ■= 


E[5i(f„,||Q‘]|-E|5i(o|Q‘]| 

E[5x„..,|[e'll 


(34) 
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Figure 8: Numerical test 1: Six realizations of the magnitude of approximate solution 
|uQg(r,X,y)| at the fixed time T = 1 with wavelength e = 1/40. In each realization, the 
central circle indicates the support of the Qol test function and the other two circles/ellipses 
indicate the supports of the initial solution. 


Here, for each wavelength, the reference solution is computed with a high level, i^et, and with 
the same Gaussian beam parameters used in all levels, i > 1. The error (34) therefore reflects 
only the stochastic collocation error, not the error in the deterministic asymptotic solver. 

Figure]^ shows the relative error, S{ri), in (34) at time T = 1, computed by the proposed 
method, versus the number of collocation points, rj, for various wavelengths. It also shows the 
convergence of the relative error in ten Monte Carlo runs, computed using the same reference 
solution as above, with rj representing the nnmber of samples. 

We observe a fast spectral convergence rate of the stochastic collocation error in the 
proposed method, due to the high stochastic regularity of the Qol. A simple linear regression 
through the data points shows that the rate of convergence of Monte Carlo is 0.45, which 
is very slow. Consequently, the decay in the stochastic collocation error is mnch faster than 
the decay in Monte Carlo error. We also note that as £ decreases, the decay rate does not 
deteriorate. This points to the existence of nniform bounds (25). 


5.2 Numerical test 2: A Lens 

In this example, the uncertain wave speed is described by a random vector y = (?/i,?/ 2 ,l/s) 
containing N = 3 independent uniformly distribnted random variables, given by 

c(x, y) = 1 — 

where 

yi ~ W(0,0.4), y 2 ~ fY(0.65, 0.85), y^ ~ 1). 
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Figure 9: Numerical test 1: Relative error Eirf) at time T — 1 versus the number of collo¬ 
cation points ri (or the number of samples in the case of Monte Carlo sampling), for various 
wavelengths. The proposed method performs a fast spectral convergence, while Monte Carlo 
sampling has a slow algebraic convergence. The rate of convergence of Monte Carlo sampling, 
obtained by linear regression through the data points, is 0.45. 
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(a) y = (0.2,0.75,0.5) 


(b) y= (0.4,0.75,0.5) 



1.5r 

1 

0.5 

^2 0 - 

-0.5 

-1 

-1.5^ 

-2 



Xi 


Xi 


(c) y = (0.2,0.75,1.0) 


(d) y= (0.4,0.75,1.0) 


Figure 10: Numerical test 2: Ray tracing solution for four realizations of y. A cusp caustic 
and two fold caustics are formed in all realizations. Contour lines of c(x, y) overlaid in red. 
The left transparent band shows where the initial amplitude is above 1/2; the right band 
shows how this band has been transported by the rays at time t = 2.5. Circle indicates the 
support of Qol test function. 


The initial data are assumed to be deterministic and given by ^ with 

$o(x) = -a;i, = 

The problem models a plane wave that is refracted by a lens of uncertain shape and strength. 

The wave speed varies in the spatial domain and caustics may consequently form. Figure 
T0| shows the ray tracing solution for four different realizations of the random vector y. A 
cusp caustic and two fold caustics are formed inside the support of the Qol test function for 
the last three realizations, but not in the first one. 

Figure [TT| shows the magnitude of the approximate solution x, y) | with wavelength 

£ = 1/20 at various time instances and a fixed random vector, y = (0.4,0.75,1), which 
corresponds to the realization in Figure [T^. The central circles indicate the support of the 
Qol test function. 
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(c) t = 2 (d) t = 2.5 

Figure 11: Numerical test 2: Magnitude of approximate solution x, y)| with wave¬ 

length £ = 1/20 at various time instances and a fixed random vector y = (0.4, 0.75,1), which 
corresponds to the realization in Figure [Toli. Circle indicates the support of Qol test function. 
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Figure 12: Numerical test 2: Relative error £(ri) at time T = 2.5 versus the number of 
collocation points, rj (or the number of samples in the case of Monte Carlo sampling), for 
various wavelengths. The proposed method performs a fast spectral convergence, while Monte 
Carlo sampling has a slow algebraic convergence. The rate of convergence of Monte Carlo 
sampling, obtained by linear regression through the data points, is 0.62. 


Figure 12 shows the relative error, £{r)), in (34) at time T — 2.5, computed by the 
proposed method, versus the number of collocation points, r} (or the number of samples in 
the case of Monte Carlo sampling), for various wavelengths. 

Similar to the first numerical test, we observe a fast spectral convergence rate of the 
stochastic collocation error in the proposed method. The convergence rate of Monte Carlo, 
given by linear regression, is in this case 0.62, and the error again decays much more slowly. 
Furthermore, the error curves of the proposed method have a rather unform shape for all the 
£ used, suggesting that (25) holds. 


6 Conclusion 

We have proposed a novel stochastic spectral asymptotic method for the forward propagation 
of uncertainty in high-frequency waves generated by highly oscillatory initial data. The source 
of uncertainty is the wave speed and/or the initial data, characterized by a finite number of 
independent random variables with known probability distributions. The proposed method 
combines a sparse stochastic collocation method for propagating uncertainty and a Gaussian 
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beam summation method for propagating high-frequency waves. Fast error convergence 
is attained only in the presence of Qols which are smooth with respect to input random 
parameters independent of the wave frequeney. 

The wave solution is highly oscillatory in both physical and stochastic spaces, and its 
derivatives clearly cannot be bounded independently of the frequency. A priori, Qols based on 
the solution would have the same behavior. However, our main result is that there are in fact 
quadratic Qols that are smooth with uniformly bounded derivatives in the stochastic space. 
Through both theoretical arguments for simplified problems and numerical experiments for 
more complicated problems, we have verified the spectral convergence of the proposed method 
for such a quadratic Qol, which represents the local wave strength. This shows that the 
proposed method may be a valid alternative to the traditional Monte Carlo method. 

Future directions include a rigorous proof of the uniform bounds in Conjecture 1 for 
the quadratic quantity considered in this work and a regularity analysis of other types of 
nonlinear Qols. 
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